Mon. Not. R. Astron. Soc. 000, 000-000 (0000) 



Printed 5 February 2008 



(MN IATeX style file v2.2) 



Dissipationless collapse of a set of N massive particles 



Fabrice Roy^)* and Jerome 

Perez (i,2)| 

M Ecole Nationale Superieure de Techniques Avancees, Unite de Mathematiques Appliquees, 32 Bd Victor, 75015 Paris, France 
( 2 ) Laboratoire de I'Univers et de ses THeories, Observatoire de Paris-Meudon, 5 place Jules Janssen, 92350 Meudon, France 



5 February 2008 



m 

O 
O 

(N 
-i— > 
O 

O 

o 

m 



oo 
O 

m 
o 

Oh 

6 

03 



ABSTRACT 

The formation of self-gravitating systems is studied by simulating the collapse of a set of 
N particles which are generated from several distribution functions. We first establish that 
the results of such simulations depend on N for small values of N. We complete a previous 
work by Aguilar & Merritt concerning the morphological segregation between spherical 
and elliptical equilibria. We find and interpret two new segregations: one concerns the 
equilibrium core size and the other the equilibrium temperature. All these features are 
used to explain some of the global properties of self-gravitating objects: origin of globular 
clusters and central black hole or shape of elliptical galaxies. 

Key words: methods: numerical, N-Body simulations - galaxies: formation - globular 
clusters: general. 



1 INTRODUCTION 

It is intuitive that the gravitational collapse of a set of N 
masses is directly related to the formation of astrophysical 
structures like globular clusters or elliptical galaxies (the pres- 
ence of gas may complicate the pure gravitational A-body 
problem for spiral galaxies) . From an analytical point of view, 
this problem is very difficult. When A is much larger than 2, 
direct approach is intractable, and since Poincare results of 
non analyticity, exact solutions may be unobtainable. In the 
context of statistical physics, the situation is more favorable 
and, in a dissipationless approximation 1 , leads to the Colli- 
sionless Boltzmann Equation (hereafter denoted by CBE) 



21 + 2L + rn 2t 2L-q 

dt ^ dr dv dp 



(1) 



where / = / (r, p, t) and V — ( r > t) are respectively the 
distribution function of the system with respect to the canon- 
ically conjugated (r, p) phase space variables and th e mean 
field gravitational potential. As noted initially by iHenorJ 
(1960), this formalism holds for such systems if and only 
if we consider A identical point masses equal to m. This 
problem splits naturally into two related parts: the time de- 
pendent regime and the stationary state. We can reasonably 
think that these two problems are not completely understood. 
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1 The dissipationless hypothesis is widely accepted in the context 
of gravitational iV-body problem because the ratio of the two-body 
relaxation time over the dynamical time is of the order of N. For 
a system composed of more than ~ 10 4 massive particles a study 
during a few hundreds dynamical times can really be considered 
as dissipationless, the unique source of dissipation being two-body 
encounters. 



The transient time dependent regim e was investigated mainly 
consid e ring self-similar solut i ons ^Ly ndon-Bell & Eselcton 
dl980l). iHenriksen fc Widrow| Jl995t) . iBlottiau et alJ dl98gl) 
and lLancellhtti^^Ktes^ffii^ ~ il200ll) ). These studies conclude 
that power law solutions can exist for the spatial dependence 
of the gravitational potential (with various powers). Never- 
theless, there is no study which indicates clearly that the 
time dependence of the solutions disappears in a few dynami- 
cal times, giving a well defined equilibriu m-like state. On the 
other hand, applying Jeans theorem (e.g. lBinnev fc Tremainel 
l|l987t) hereafter BT87, p. 220), it is quite easy to find a sta- 
tionary solution. For example, every positive and integrable 
function of the mean field energy per mass unit E is a potential 
equilibrium distribution function for a spherical isotropic sys- 
tem. Several approaches are possible to choose the equilibrium 
distributio n function. Thermo d ynamics (Violen t Relaxation 
parad igm: iLvnden-Beul jl967lh IChavanisI l)2002lh iNakamural 
(2000)) indicate that isothermal spheres or polytropic systems 
are good candidates. Stability analys es can be spl i t into two 
categories. In the CBE context (see IPerez fc Alvl l)l998h for 
a review), it is well known that spherical systems (with de- 
creasing spatial density) are generally stable except in the case 
where a large radial anisotropy is present in the velocity space. 
Thi s is the Radia l Orbi t Inst ability, hereafter de noted by ROI 
(see lPerez fc Alvl dl998l) . and lPerez et alJ dl998l) for a detailed 
analytic and numeric study of these phenomena) which leads 
to a bar-like equilibrium state in a few dynamical times. In 
the context of thermo dynamics of self -gravitating systems, in 
a pioneering work bv lAntonovl l|l962lh it was shown that an 
important density contrast l eads to the collapse of the core of 
system (see IChavanisI (|^)03) for details). 

In all these studies there is no definitive conclusion, and the 
choice of the equilibrium distribution remains unclear. Intro- 
ducing observations and taking into account analytical con- 
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straints, several models are possible: chronologically, we can 
cite (see for example BT87, p. 223-239) the Phimmer model 
(or other polytropic models), de Vaucouleurs r 1//4 law, King 
and isochrone Henon model or more recently the very sim- 
ple but interesting Hernquist model l|Hernauistl l)l990t0 for 
spherical isotropic systems. In the anisotropic case, Ossipkov- 
Merritt or generalized polytropes can be considered. Finally 
for non spherical systems, there also exists some models re- 
viewed in BT87 (p. 245-266). Considering this wide vari- 
ety of possibilities, one can try to make accurate numeri- 
cal simulations to clarify the situation. Surprisingly, such a 
program has not been comp letely carried on. In a pioneer- 
ing work, Ivan Albadal il982T) remarked that the dissipation- 
less collapse of a clumpy cloud of iV equal masses could 
lead to a final stationary state that is quite similar to el- 
liptical galaxies. T his kind of study was rec onsidered in an 
important work bv lAeuilar fc Merrittl (^.990), with more de- 
tails and a crucial remark concerning the correlation between 
the final shape (spherical or oblate) and the virial ratio of 
the initial state. These authors explain this feature invokin g 
ROI. Some more r ecen t studies llCan nizzo fc Hollisteil 1^)92) , 
iBoilv et al.l Jl999ft and lTheis fc Spurzeml \\99$) ) concentrate 
on so me particularit i es of the preceding works. Fina l ly, tw o 
works JPantas et alJ J20021) and lCarpintero fc Muzziol JlQ95l> s ) 
develop new ideas considering the influence of the Hubble flow 
on the collapse. However, the problem is only partially de- 
picted. 

The aim of this paper is to analyse the dissipationless collapse 
of a large set of N Body systems with a very wide variety of 
'realistic' initial conditions. As we will see, the small number 
of particles involved, the numerical technique or the specificity 
of the previous works did not allow their authors to reach a 
sufficiently precise conclusion. The layout of this paper is as 
follows. In section 2 we describe in detail the numerical proce- 
dures used in our experiments. Section 3 describes the results 
we have obtained. These results are then interpreted in section 
4, where some conclusions and perspectives are also proposed. 



2 NUMERICAL PROCEDURES 
2.1 Dynamics 

The Treecode used to perform our simulations is a modified 
version of the lBarnes fc Hull 1^986) Treecode, parallelised by 
D. Pfenniger using the MPI library. We implemented some 
computations of observables and adapted the code to suit our 
specific problems. The main features of this code are a hi- 
erarchical 0(N \og(N)) algorithm for the gravitational force 
calculation and a leap-frog algorithm for the time integration. 
We introduced an adaptative time step, based on a very simple 
physical consideration. The time step is equal to a fraction nts 
of the instantaneous dynamical time Td ( 2 ) , i.e. At = Td/nts 



2 The fraction nts is adapted to the virial parameter n and ranges 
roughly from nts = 300 when n = 90 to nts = 5000 when n = 08. 
The dynamical time we used is given by 



i=l 

N 



. The simulations were run on a Beowulf cluster (25 dual CPU 
processors whose speed ranges from 400MHz to 1GHz). 



2.2 Initial Conditions 

The initial virial ratio is an important parameter in our sim- 
ulations. The following method was adopted to set the virial 
ratio to the value VinMai- Positions r,: and velocities v< are 
generated. We can then compute 

2K 



U 



where 



1 

K = ^2 2 miV * 2 



and 



U = 



N 



2 ~^ {max^n-Tj) . 



2 6 2))V2 



(2) 



(3) 



(4) 



In this relation e is a softening parameter whose value is dis- 
cussed in section T2. 3. 21 As the potential energy depends only 
on the positions, we obtain a system with a virial ratio equal 
to Vi n itiai just by multiplying all the particle velocities by the 
factor {Vinitiai /Vp) 1 ^ 2 . For convenience we define 



V = | Vinitiai | x 10 2 



2.2.1 Homogeneous density distribution (H v ) 



(•») 



As we study large A-body systems, we can produce a ho- 
mogeneous density by generating positions randomly. These 
systems are also isotropic. We produce the isotropic velocity 
distribution by generating velocities randomly. 



2.2.2 Clumpy density distribution (C^) 

A type of inhomogeneous systems is made of systems with a 
clumpy density distribution. We first generate n small homoge- 
neous spherical systems with radius R g . Centers of these sub- 
systems are uniformly distributed in the system. The empty 
space is then filled using a homogeneous density distribution. 
In the initial state, each clump contains about 1% of the total 
mass of the system and have a radius which represents 5% 
of the initial radius of the whole system. These systems are 
isotropic. 



2.2.3 Power law r a density distribution (P%) 

We first generate the f> and z cylindrical coordinates using 
two uniform random numbers ui and u^: 

[z, <p) = (2ui - 1, 27TU 2 ) . (6) 

Using the inverse transformation method, if 



r = RF' 1 (u) with F (r) = ^ I x 2 ' a dx 



(7) 



where R is the radius of the system, it is a uniform random 
number, i< 1 and 



S = I x 2 a dx 



(8) 
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then the probability density of r is proportional to r 2 a , and 
the mass density p is proportional to r~ a . Finally, one gets 



z z cos If 



- z A sm ip 
vz 



(9) 



These systems are isotropic. 



2.2.4 Gaussian velocity distribution(G^ ) 

Most of the systems we use have a uniform velocity distribu- 
tion. But we have also performed simulations with systems 
presenting a gaussian initial velocity distribution. These sys- 
tems are isotropic, but the x-, y- and z-components of the 
velocity are generated following a gaussian distribution. Using 
a standard method we generate two uniform random numbers 
ui and U2, and set 



Vi — y— 2a 2 In ui cos(27rcr 2 it2) 



i = x,y,z 



(10) 



where a is the gaussian standard deviation. 
2.2.5 Global rotation (R s v ) 

Some of our initial systems are homogeneous systems with a 
global rotation around the Z-axis. The method we choose to 
generate such initial conditions is the following. We create a 
homogeneous and isotropic system (an H-type system). We 
then compute the average velocity of the particles. 



N 

j=i 



(11) 



We project the velocities on a spherical referential, and add a 
fraction of v to v$ with regard to the position of the particle. 
We set 



, f P iV 



(12) 



where / is a parameter of the initial condition, pi is the dis- 
tance from the particle to the Z-axis and R the radius of the 
system. The amount of rotation induced by this method can 
be evaluated through the ratio: 



H = K ro t/K, 



(13) 



where K is the total kinetic energy defined above, 
whereas K ro t is the r otation kinetic energy defined by 
iNavarro fc White I dl993l) : 



r , 1 v~* (J* ' Jtot) 2 

2 tT [r 2 -( r! .J tot ) 2 ] 



(14) 



Above, Ji is the specific angular momentum of particle i, Jtot 
is a unit vector in the direction of the total angular momentum 
of the system. In order to exclude counter-rotating particles, 
the sum in equation 114H is actually carried out only over those 
particles which verify the condition (Ji ■ Jtot) > 0. 



The number of particles whose mass is M ^ m ^ M + dM 
is n(M)dM. In some models, the value of a and j3 depends 
on the range of mass that is considered. We have used several 
types of mass functions, among them the initial mass function 
given by iKrouoal (|200lT ) (fc = I), the one given bv ISaloeteil 
(1955) (fc = II) and an M _1 mass function (k = III). In 
order to generate masses following these functions, we first 
calculate at to produce a continuous function. We then can 
calculate the number of particles whose mass is between M 
and M + dM. We generate n(M) masses 



rrn = M + udM 



1 < i < n(M) 



(16) 



where ^ u ^ 1 is a uniform random number. In the initial 
state, these systems have a homogeneous number density, a 
quasi homogeneous mass density and they are isotropic. 

2. 2. 7 Nomenclature 

We indicate below the whole set of our non rotating initial 
conditions. 

• Homogeneous models: H 8 8, H 79 , H 6 o, H 50 , H 4 o, H 30 , 
H20, H15 and H10 

• Clumpy C^models: C|?, C§g, C|?, Cg, C§g, Cgg, Cf A , 
C10, C07 and C10 

. Power Law models: P 2 °, P 2 9 °, P^ °, P° 5 , P 10 °, Pj 8 5 
and Pi 5 

• Gaussian velocity profiles models: G50, G 2 q, G50, G12 
and Gi 

• Mass spectra M* models: M' 50 , M50, M5", M35, M25, 
Mli 1 and M^ 7 

For all these models we ran the numerical simulations 
with 30 000 particles (see § lO 



2.3 Observables 

2.3.1 Units 

Our units are not t he commonly used ones (see 
iHeggie fc Mathieul lll98d) '). We did not set the total en- 
ergy E of the system to —0.25 because we wanted to prescribe 
instead the initial virial ratio Visual, the size 7? of the system 
and its mass M. We thus have M = 1, R = 10 and G = 1, 
and values of VinUial and E depending on the simulation. We 
can link the units we have used with more standard ones. We 
have chosen the following relationships between our units of 
length and mass and common astrophysical ones: 



M = 10" M Q and R = 10 pc 
Our unit of time ut is given by: 



"' ^ lit- 4 - 72 1011 S = L5 ° * 



(17) 



(18) 



where variables X s are expressed in our simulation units and 
variables X c in standard units. 



2.2.6 Power-law initial mass function (Ml) 

Almost all the simulations we made assume particles with 
equal masses. However, we have created some initial systems 
with a power-law mass function, like 

n(M) = aM 13 (15) 



2.3.2 Potential softening and energy conservation 

The non conservation of the energy during the numerical evo- 
lution has three main sources. 

The softening parameter e introduced in the potential calculus 
(cf. equation is an obvious one. This parameter introduces 
a lower cutoff A in the resolution of length in the simulations. 
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Following iBarnes fc Hull dl989l) . structural details up to scale 
A < lOe are sensitive to the value of e. Moreover, in order to 
be compatible with the collisionless hypothesis, the softening 
parameter must be greater than the scale where impo rtant 
collisions can occur. Still following IBarnes fc Hud l)l98fll) . this 



10 



G{m) 
(v 2 ) 



(19) 



In our collapse simulation with 3 • 10 4 particles, this results 
in £77 > 2/3. The discretization of time integration introduces 
inevitably another source of energy non conservation, particu- 
larly during the collapse. The force computation also generates 
errors. The choice of the opening angle O, which governs the 
accuracy of the force calculation of the Treecode, is a compro- 
mise between speed and accuracy. For all these reasons, wc 
have adopted e = 0.1. This choice imposes r/ > 6 (for 30 • 10 4 
particles). This trade-off allowed to perform simulations with 
less than 1% energy variation without requiring too much com- 
puting time. For each of our experiments, the total CPU time 
ranges between 3 to 24 hours for 3000 u t and 3 • 10 4 particles. 
The total agregated CPU time of all our collapse experiments 
is approximately 6 months. 

We have tested two other values of the softening parameter 
(e = 0.03 and e = 0.3) for several typical simulations. These 
tests did not reveal significant variations of the computed ob- 
servables. 



2.3.3 Spatial indicators 

As indicators of the geometry of the system, we computed 
axial ratios, radii containing 10% (-Rio), 50% (Rbo) and 90% 
(-R90) of the mass, density profile p(r) and equilibrium core 
radius. The axial ratios are computed with the eigenvalues Ai, 
A2 and A3 of the (3x3) inertia matrix J, where A3 ^ A2 ^ Ai 
and, if the position of the particle i is = (xij ; X2,% ; 2:3,1) 



JV 



for p / v = 1, 2, 3 



for p = 1, 2, 3 



(20) 



The axial ratios a\ and 02 are given by 02 = A1/A2 1.0 and 
ai = A3/A2 < 1.0. 

The density profile p, which depends only on the radius r, 
together with the Rs (8 — 10, 50, 90) have a physical mean- 
ing only for spherical or nearly spherical systems. For all the 
spatial indicators computations we have only considered par- 
ticles whose distance to the center of mass of the system is less 
than 6 x R50 of the system. This assumption excludes particles 
which are inevitably 'ejected' during the collapse 3 . 
After the collapse a core-halo structure forms in the system. 
In order to measure the radius of the core, we have com- 
put ed the density rad i us as defined by Casertano and Hut 
fsee ICasertano fc Hutl l)l985h '). The density radius is a good 
estimator of the theoretical and observational core radius. 
We have also computed the radial density of the system. The 



3 The number of excluded particles ranges from 0% to 30% of the 
total number of particles, depending mostly on n. For example, the 
number of excluded particles is 0% for Hgo, 3% for C|^, 5% for Hso, 
22% for C 2 ° and 31% for H i0 . 



density is computed by dividing the system into spherical bins 
and by calculating the total mass in each bin. 

2.3.4 Statistical indicators 

When the system has reached an equilibrium state, we com- 
pute the temperature of the system 



T 



2(K) 



(21) 



3Nk B 

where K is the kinetic energy of the system, ks is the Boltz- 
mann constant (which is set to 1) and the notation (A) denotes 
the mean value of the observable A, defined by 

N 

i=i 

In order to characterise the system in the velocity space we 
have computed the function 



K (r) = 



2 ( V lrad) r<r . <r _i 
( V i,tari) 



(23) 



r ^ r i <r+dr 



where Vi, r ad is the radial velocity of the i particle, and Vijan 
its tangential velocity. For spherical and isotropic systems 
(ai ~ 02 — 1 and n (r) ~ 1), we have fitted the density by 
1- a polytropic law 



p = palp 1 

which corresponds to a distribution function 

/ (E) oc E"<-' i/2 

2- an isothermal sphere law 

lb Is 2 

p = pie v/ 

which corresponds to a distribution function 



f(E) = 



Pi 



(2%s 



2^3/2 



(24) 

(25) 
(26) 
(27) 



Using the least square method in the ln(p) — ln(tp) plane we 
get (7,ln(p )) and (s 2 , In (pi)) . 



3 DESCRIPTION OF THE RESULTS 

We have only studied systems with an initial virial ratio cor- 
responding to 7] € [7,88]. In such systems, the initial velocity 
dispersion cannot balance the gravitational field. This pro- 
duces a collapse. After this collapse, in all our simulations the 
system reaches an equilibrium state characterised by a tempo- 
ral mean value of the virial ratio equal to —1, i.e 77 = 100, and 
stationary physical observables. These quantities (defined in 
section 2) are presented in a table of results in the appendix 
of this paper. The following results will be discussed and in- 
terpreted in section [1] 



3.1 Relevant number of particles 

In all previous w o rks on this subject ifva n Alba dal <ll982|) . 
Aguilar fc Merrittl (|l99Ctl . ICannizzo fc Hollisterl ljl992h and 
Bolh^t^Lri^99li '). the authors did not really consider the 
influence of the number of particles on their results. In the 
first two and more general works, this number is rather small 
(not more than a few thousands in the largest simulations). 
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Figure 2. Axial ratios of equilibrium states reached from Homoge- 
neous, Clumpy, Gaussian velocity dispersion, Power law and Mass 
spectrum initial conditions. 



of particles used in each case ranges from 10 2 to 10 5 . We can 
see in Figure Q that some observables are N-dependent when 
N < 3.10 4 . In particular, R50 and Rd present a monotonic 
variation larger than 50% when N varies from 10 2 to 10 s and 
the ellipticity of the final state is overestimated for small val- 
ues of N. As a conclusion of this preliminary study, we claim 
that the relevant number of particles for collapse simulations is 
N ^ 3.10 4 . As all simulations have been completed with a to- 
tal energy loss smaller than 1%, we state moreover that this re- 
sult is independent of the numerical scheme used (Treecode or 
Direct N-Body). As a consequence, the simulations presented 
hereafter have been performed using N — 3.10 4 particles. 



4 r 0.0, 



Figure 1. Influence of the number of particles on the physical ob- 
servables of a collapsing system. Axial ratios are on the top panel, 
radius containing 50% of the total mass and density radius are on 
the middle pannel and the best s 2 and 7 fit for respectively isother- 
mal and polytropic distribution function are on the bottom panel. 
All cases are initially homogeneous with 77 = 10 (solid line), 77 = 50 
(dotted line) and 77 = 80 (dashed line). The number N of particles 
used is in units of 10 3 . 

The two other studies are more specific and use typically 10 4 
and, in a few reference cases, 2.10 4 particles. In order to test 
the influence of the number of particles on the final results, we 
have computed several physical observables of some collapsing 
systems with various numbers of particles. The results are pre- 
sented in Figure Q In order to check the influence of N in the 
whole phase space, we have studied positions and velocities 
related observables: 01, 02, Rio, R50 and Rgo and parame- 
ters of isothermal and polytropic fit models namely 7 and s 2 . 
Moreover, in order to be model-independent, we have studied 
three representative initial conditions: Hso, H50 and H10, i.e. 
respectively initially hot, warm and cold systems. The number 



3.2 Morphological segregation 

An important study bv lAguilar fc Merritd (Il990l) shows that, 
in the case of an initial density profile p oc r _1 , the shape of 
the virialized state depends on 77: a very small 77 leads to a 
flattened equilibrium state, when a more quiet collapse pro- 
duces a spherical one. Our investigations concern a wide range 
of different initial conditions (homogeneous, clumpy,...) and 
show that the influence of 77 depends on the properties of the 
initial system 4 . Figure |5] shows the axial ratios of the equi- 
librium state reached by our simulations. In fact, only a few 
simulations produced a final state with an ellipticity greater 
than El. Every homogeneous initial condition (i.e. H,,, and 
resulted in a spherical equilibrium state independently of 
the values of 77 we tested. Cold clumpy systems have a weakly 
flattened equilibrium state. The only final systems with an el- 
lipticity significantly greater than El are those produced by 
the collapse of cold P™. 

Previous studies invoked ROI to explain this morphological 
segregation. However, it seems that ROI requires inhomo- 
geneities near the center to be triggered. 



4 The particular case of rotating initial conditions is discussed in a 
special section. 
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Figure 3. Density profile for models plotted in units of R50 



Figure 5. Density profile for models plotted in units of Rso- 





Figure 4. Density profile for C™ models plotted in units of -R50. 
The dashed line corresponds to the C^q model 



Figure 6. Density profile for models plotted in units of -Rso- 



3.3 Characteristic size segregation 

In addition to the morphological segregation, presented in the 
previous section, we discovered a finer phenomenon. 
In the Figures 0] and we have plotted the mass 
density of all equilibrium states produced by the collapse of 
our initial conditions as a function of the ratio r/Rso. These 
plots represent the density at the end of our simulations (after 
about 100 crossing times) . These functions do not significantly 
evolve after the collapse except for Mq 7 . For this special case, 
a comparative plot is the subject of Figure |H1 
All equilibrium states we obtain clearly fall into two cate- 
gories: 

• Flat Core Systems 
All these systems present a core halo structure, i.e. a large 
central region with a constant density and a steep envelope. 



These systems are typically such that In (Rso/Rd) < 0.5 and 
\n(R w /Rd) < -0.05. 

• Small Core Systems 
For such systems, the central density is two order of magnitude 
larger than for Flat Core systems. There is no central plateau 
and the density falls down regularly outward. These systems 
are typically such that \n(Rso/Rd) > 0.7 and ln(7?io/i?d) > 
0.1. 

The diagram In (Rio /Rd) vs ln(i?5o/-Rd) is the subject 
of Figure |§] One can see in this figure that each equilibrium 
state belongs to one or the other family except in a few par- 
ticular cases. In the Flat Core family we found all H,,, and 
Mjj systems except M7, and two P" systems namely P50 and 
P50. These systems are all initially homogeneous or slightly 
inhomogeneous (e.g. P55 5 and P^q systems). In the Small Core 
family, we found all C^° and the P§5° and P09 systems. These 
systems are all initially rather very inhomogeneous. Finally, 
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Figure 7. Density profile for M* models plotted in units of Rsq. 



Figure 9. The core size segregation: ln(flio/r c ) vs In (R50 A*e) is 
plotted for all non rotating systems. 
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Figure 8. Comparison between the evolution of the mass density 
with respect to time for C 2 q (top-panel) and Mq 7 . For each case, 
plotted times are 10,20,30,40,50,75 and 100 T d 



there are 5 systems in-between the two categories: C10, P45 5 , 
Pog 5 , P10 and Mq 7 . This last model is the only one which mi- 
grates from Flat core set (when t ~ lOTd) to the edge of the 
small core region (when t ~ lOOTd). 



3.4 Equilibrium Distribution Function 

In order to compare systems in the whole phase space, we 
fitted the equilibrium state reached by each system with two 
distinct isotropic models, e.g. polytropic and isothermal (see 
equations or BT87, p. 223-232). Figure HT71 shows these 
two fits for the P55 5 simulation. The technique used for the 
fit is described in section 2 of this paper. The result obtained 
for this special study is the following: the equilibrium states 
reached by our initial conditions can be fitted by the two 
models with a good level of accuracy. As long as f) < 70, 
the polytropic fit gives a mean value 7 = 4.77 with a stan- 
dard deviation of 2.48 1CP 1 . This deviation represents 5.1% 
of the mean value. This value corresponds typically to the well 
known Plummer model for which 7 = 5 (see BT87 P. 224 for 
details). When the collapse is very quiet ( typically r\ > 70 ) 
polytropic fit is always very good but the value of the index is 
much larger than Plummer model, e.g. 7 = 6.86 for H79 and 
7 = 7.37 for Hg8. The corresponding plot is the subject of the 
Figure [T2"l All the data can be found in appendix. As we can 
see on the example plotted in Figure E3 the isothermal fit is 
generally not as good as the polytropic one. On the whole set 
of equilibria, isothermal fits give a mean value s 2 = 2.5 10" 2 
with a standard deviation of 1.6 10 -2 (60%). The correspond- 
ing plot is the subject of Figure All the data can be found 
in appendix. 

In fact both isothermal and polytropic fits are reasonable: 
as long as the model is able to reproduce a core halo structure 
the fit is correct. The success of the Plummer model, which 
density is given by 

„-, -5/2 



p{r) 



4Trb 3 



1 + 



can be explained by its ability to fit a wide range of models 
with various ratio of the core size over the half-mass radius. 
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Figure 10. Polytropic and isothermal fit for the Pgj} 5 simulation. 
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Figure 11. Best fit of the s 2 parameter for an isothermal model 
for all non rotating systems studied. The error bar correspond to 
the least square difference between the data and the model. 




Figure 12. Best fit of the 7 parameter for a polytropic model for 
all non rotating systems studied. The error bars correspond to the 
least square difference between the data and the model. 
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The adjustment of this ratio is made possible by varying the 
free parameter b. We expect that other core halo models like 
King or Hernquist models work as well as the Plummer model. 
As a conclusion of this section, let us say that as predicted 
by theory there is not a single universal model to describe the 
equilibrium state of isotropic spherical self-gravitating system. 

3.5 Influence of rotation 

We saw in section l3~!?l a source of flattening for self-gravitating 
equilibrium. Let us now show the influence of initial rotation, 
which is a natural candidate to produce flattening. The way 
we have added a global rotation and the significance of our 
rotation parameters / and \x are explained in section 12.2.51 
The set of simulations performed for this study contains 31 
different elements. The initial virial ratio ranges from r\ — 
10 to 77 = 50, and the rotation parameter from f — (i.e. 
fi = 0) to / = 20 (i.e. /i = 0.16 when 77 = 50). As a matter 
of fact, equilibrium states always preserve a rather important 



Figure 13. Axial ratios as for different values of r\ as a function of 
the initial solid rotation parameter / 



part of the initial rotation 5 and, observed elliptical gravitating 
syst ems generally posses s very small amount of rotation (see 
e.g. ICombes et all l)l995h 'l. We thus exclude large values of / 

Our experiments exhibit two main features (see Figure tTsT l: on 
the one hand, rotation produces a flattened equilibrium state 
only when / exceeds a triggering value (typically / = f ~ 4) . 
On the other hand, we have found that for a given value of 
77, the flatness of the equilibrium is roughly /—independent, 
provided that / > f a . 



5 We observed that /i is always smaller in the equilibrium state than 
in the initial one, typically each rotating systems conserves 65% of 
the initial /i 
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Figure 14. Evolution of the temperature as a function of time 
3.6 Thermodynamical segregation 

As we study isolated systems, the total energy E contained in 
the system is constant during the considered dynamical evo- 
lution . This property remains true as long as we consider 
collisionless evolutions. For gravitational systems, this means 
that we cannot carry out any simulation of duration larger 
than a few hundred dynamical times. We have obviously taken 
this constraint into account in our experiments. All systems 
which experience a violent relaxation reach an equilibrium 
state which is stationary in the whole phase space. Spatial be- 
haviour like morphological segregation produced by ROI was 
confirmed and further detailled thanks to our study. A new 
size segregation was found in section HOI Now let us consider 
another new segregation appearing in the velocity space. Each 
equilibrium state is associated with a constant temperature T, 
calculated using equation 12111 . More precisely, we have calcu- 
lated the temporal mean value 6 of the temperatures, evalu- 
ated every one hundred time units. As we can see in Figure 
1141 after the collapse and whatever the nature of the initial 
system is, the temperature is a very stable parameter. 
Figure Il5l shows the E — T diagram of the set of all non rotat- 
ing simulations. It reveals a very interesting feature of post- 
collapse self-gravitating systems. 

On the one hand, the set of systems with a total energy 
E > —0.054 can be linearly fitted in the E — T plane. We call 
this set Low Branch 1 (hereafter denoted by LB1, see Figure 
I15JI . On the other hand, the set of the systems with a total 
energy E < —0.054 splits into two families. The first is an 
exact continuation of LB1. Hence we named it Low Branch 2 
(hereafter denoted by LB2). The second can also be linearly 
fitted, but with a much greater slope (one order of magnitude) . 
We label this family High Branch (hereafter denoted by HB). 
In LB1 or LB2, we find every H, G, and M systems with 
77 > 25, every P and every C with n > 10. In HB, we find C?o 
and every H, G and M systems with r\ < 25. This segregation 
thus affects violent collapses (cold initial data): on the one 
hand, when 77 > 25 all systems are on LB1, on the other hand 



6 The temporal mean value is computed from the time when the 
equilibrium is reached until the end of the simulation 
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Figure 15. Energy- Temperature diagram 

for r\ < 25, initially homogeneous or quasi- homogeneous (e.g. 
C?q) systems reach HB when initially inhomogeneous systems 
stay on LB2 instead. 



4 INTERPRETATIONS, CONCLUSIONS AND 
PERSPECTIVES 

Let us now recapitulate the results we have obtained and pro- 
pose an interpretation: 

(i) The equilibrium state produced by the collapse of a set 
of N gravitating particles is N— independent provided that 
N > 3.0 10 4 . 

(ii) Without any rotation, the dissipationless collapse of a 
set of gravitating particles can produce two relatively distinct 
equilibrium states: 

• If the initial set is homogeneous, the equilibrium has a 
large core and a steep envelope. 

• If the initial set contains significant inhomogeneities 
(n > 10 for clumpy systems or a > 1 for power law sys- 
tems) , the equilibrium state has only a small core around of 
which the density falls down regularly. 

The explanation of this core size segregation is clear: it is 
associated to the Antonov core-collapse instability occurring 
when the density contrast between central and outward re- 
gion of a gravitating system is very big. As a matter of fact, 
if the initial set contains inhomogeneities, these collapse much 
more quickly than the whole system 7 and fall quickly into the 
central regions. The density contrast becomes then very large 
and the Antonov instability triggers producing a core collapse 
phenomenon. The rest of the system then smoothly collapses 
around this collapsed core. If there are no inhomogeneities in 
the initial set, the system collapses as a whole, central den- 
sity grows slowly without reaching the triggering value of the 
Antonov instability. A large core then forms. Later evolution 
can also produce core collapse: this is what occurs for our Mq 7 
system (see Figurc0 . This is an initially homogeneous system 

7 Because their Jeans length is much more smaller than the one of 
the whole system. 
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with Kroupa mass spectrum which suffers a very strong col- 
lapse. As the mass spectrum is not sufficient to bring quickly 
enough a lot of mass in the center of the system, Antonov 
instability does not trigger and a large core forms. As the col- 
lapse is very violent, an increasing significant part of particles 
are progressively "ejected" and the core collapse takes pro- 
gressively place. This is the same phenomenon which is gener- 
ally invoked to explain the collapsed c ore of some old globular 
clusters fe.g. iDiorgowski et al.l Tl986h '): during its dynamical 
evolution in the galaxy, some stars are tidally extracted from 
a globular cluster, to compensate this loss the cluster concen- 
trate its core, increasing then the density contrast, triggering 
sooner or later the Antonov instability. 

(iii) Without any rotation, the collapse (violent or quiet) 
of an homogeneous set of gravitating particles produces an 
E0 (i.e. spherical) isotropic equilibrium state. There are two 
possible ways to obtain a flattened equilibrium: 

• Introduce a large amount of inhomogeneity near the 
center in the initial state, and make a violent collapse (7/ < 
25). 

• Introduce a sufficient amount (/ > 4) of rotation in the 
initial state. 

These two ways have not the same origin and do not produce 
the same equilibrium state. 

In the first case, one can reasonably invoke the Radial Orbit 
Instability: as a matter of fac t, as it is explained in a lot of 
works (see IPerez et al.l l)l998l) for example) two features are 
associated to this phenomenon. First of all, it is an instability 
which needs an equilibrium state from which it grows. Sec- 
ondly, it triggers only when a sufficient amount of radial or- 
bits are present. The only non rotating flattened systems we 
observed just combine this two conditions: sufficient amount 
of radial orbits because the collapse is violent and something 
from which ROf can grow because we have seen in the previous 
point that inhomogeneities collapse first and quickly join the 
center. The fact that cold systems are more flattened than 
ones is in complete accordance with our interpretation: 
as a matter of fact, by construction, power law systems have 
an initial central overdensity, whereas clumpy systems create 
(quickly but not instantaneously) this overdensity bringing the 
collapsed clumps near the center. The ROI flattening is oblate 
(d2 ~ I and ai < 1). 

The rotational flattening is more natural and occurs when 
the centrifugal force overcomes the gravitational pressure. The 
rotational flattening is prolate ( 02 > 1 and a\ ~ I). We no- 
tice that initial rotation must be invoked with parsimony to 
explain the ellipticity of some globular clusters or elliptical 
galaxies. As a matter of fact, these objects are very weakly 
rotating systems and our study has shown that the amount of 
rotation is almost constant during the collapse. 

(iv) Spherical equilibria can be suitably fitted by both 
isothermal and polytropic laws with various indexes. It sug- 
gests that any distribution function of the energy exhibiting 
an adaptable core halo structure (Polytrope, Isothermal, King, 
Hernquist,...) can suitably fit the equilibrium produced by the 
collapse of our initial conditions. 

(v) There exists a temperature segregation between equilib- 
rium states. It concerns only initially cold systems (i.e. systems 
which will suffer a violent collapse): for such systems when r\ 
decreases, the equilibrium temperature T increases much more 
for initially homogeneous systems than for initially inhomogc- 
neous systems. On the other hand, whatever their initial ho- 



mogeneity, quiet collapses are rather all equivalent from the 
point of view of the equilibrium temperature: T increases in 
the same way for all systems as rj decreases. This feature may 
be the result of the larger influence of the dynamical friction 
induced by the primordial core on the rapid particles in a vi- 
olent collapse. 

All these properties may be directly confronted t o phys- 
ical d ata from globular clusters (see Harris catalogue lHarrisI 
(1996)) or galaxies observations. 

As a matter of fact, in the standard "bottom-up" scenario of 
the hierarchical growth of structures, galaxies naturally form 
from very inhomogeneous medium. Our study then suggests 
for the equilibrium state of such objects a potential flattening 
and a collapsed core. This is in very good accordance with the 
EO to E7 observed flatness of elliptical galaxies and may be a 
good explanation for t he presence of massiv e black hole in the 
center of galaxies fsee ISchodel et alJ (|2002)). 
On the other hand Globular Clusters observations show that 
these are spherical objects (the few low flattened clusters all 
possess a low amount of rotation), and that their core is gen- 
erally not collapsed (the collapsed core of almost 10% of the 
galactic Globular Clusters can be explained by their dynami- 
cal evolution through the galaxy). Our study then expect that 
Globular Clusters form from homogeneous media. 
These conclusions can be tested using the E — T plane. As 
a matter of fact, we expect that an E — T plane build from 
galactic data would not present any High Branch whereas the 
same plane build from Globular Clusters data would. 
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Table Al. Homogeneous Initial Conditions: H, 
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Table A4. Mass Spectrum Initial Conditions: 
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0.4 


99 


1.01 


0.98 


1.18 


2.51 


225.5 


1.54 


4.66 


-2 


3.59 


-3 


4.08 


5.23 


35 


Krou 


0.2 


98 


1.01 


0.99 


1.55 


3.18 


39.8 


2.09 


4.66 


-2 


2.97 


-3 


3.51 


4.95 


51 


1/M 


0.2 


95 


1.01 


0.98 


1.79 


4.19 


14.2 


2.83 


4.67 


-9 


2.39 


-4 


3.15 


4.48 


50 


Krou 


0.1 


96 


1.02 


0.98 


1.93 


4.09 


15.1 


2.87 


4.60 


-1 


2.45 


-3 


3.19 


4.50 


50 


Salp 


0.1 


96 


1.01 


0.98 


2.03 


4.13 


15.1 


2.87 


4.70 


-2 


2.45 


-2 


3.08 


4.49 



Table A5. Gaussian Velocity Dispersion Initial Conditions: G^ 



V 


a 


A E 


V 


Ol 


a.2 


Rio 


Rso 


R90 


Rd 


7 


E 2 


s 2 


S 2 2 


T 


—E 
























(X 10 2 ) 


(X 10 2 ) 




(X 10 ) 


(X 10 11 ) 


(X 10 2 ) 


48 


Gl 


0.0 


95 


1.00 


1.00 


1.72 


3.98 


14.6 


2.23 


4.66 


-5 


2.64 


-3 


3.13 


4.50 


49 


G2 


0.0 


95 


1.02 


0.99 


1.74 


4.02 


14.4 


2.28 


4.65 


-6 


2.61 


-3 


3.13 


4.50 


50 


G3 


0.0 


95 


1.00 


1.00 


1.90 


4.08 


14.1 


2.56 


4.72 


-8 


2.52 


-2 


3.09 


4.50 


12 


G4 


0.4 


118 


1.03 


0.98 


0.56 


1.77 


1312.0 


0.64 


4.56 


-1 


5.33 


-10 


6.08 


5.63 


50 


G5 


0.0 


96 


1.01 


0.99 


1.98 


4.11 


14.8 


2.72 


4.71 


-1 


2.50 


-2 


3.19 


4.50 
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